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Abstract: Inteipretation of multi-angle spectro-polarimetric data in remote sensing of atmospheric 
aerosols require fast and accurate methods of solving the vector radiative transfer equation (VRTE). The 
single and double scattering approximations could provide an analytical framework for the inversion 
algorithms and are relatively fast, however accuracy assessments of these approximations for the aerosol 
atmospheres in the atmospheric window channels have been missing. This paper provides such analysis 
for a vertically homogeneous aerosol atmosphere with weak and strong asymmetry of scattering. In both 
cases, the double scattering approximation gives a high accuracy result (relative error ~ 0.2%) only for 
the low optical path ~ 10”^. As the error rapidly grows with optical thickness, a full VRTE solution is 
required for the practical remote sensing analysis. It is shown that the scattering anisotropy is not 
important at low optical thicknesses neither for reflected nor for transmitted polarization components of 
radiation. 
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1. Introduction 

A complete solution of the vector radiative transfer equation (VRTE) is time-consuming. At the 
same time, interpretation of the airborne (Chowdhary et al., 2001) as well as of the cun'ent (Dubovik et 
a/., 201 1) and future satellite multi-angle spectro-polarimetric measurements requires fast and sufficiently 
accurate methods of VRTE solution. The polarimetric accuracy of the POLDER-type instrument is -2% 
(http://smsc.cnes.fr /POLDER/ A_calibration.htm). The more precise systems target the accuracy of 
-0.2% for "the measured orthogonal polarization states" (Mishchenko et al., 2007) or the degree of linear 
polarization (Diner et al. 2007, 2010). 

A fast approximate VRTE solution can be obtained using the first two orders of scattering. These 
approximations offer convenience by providing an analytical framework for the inversion algorithms. The 
original derivation for the first two scattering orders for Rayleigh scattering was given by Hammad and 
Chapman (1939), Hammad (1948) with numerical results presented in (Hammad, 1953). Later, Dave 
(1964) gave a solution and numerical analysis for the first three orders of Rayleigh scattering. In his 
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study, Dave (1964, p.313) concluded that it is sufficient to consider only primary and secondary 
scattering when the scattering optical thickness is 0.05. But the number of necessary orders of scattering 
increases very rapidly with optical thickness and is about three and eight for optical thickness 0.1 and 
1.0, respectively" or higher to achieve 0.1% level of accuracy. 

Using the adding-doubling method, Hovenier (1971) derived equations for the first few orders of 
scattering for the plane-parallel homogeneous atmosphere. This work used an analytical integration over 
the optical thickness and numerical integration over the angle. From analysis for the Rayleigh and two 
types of Henyey-Greenstein phase matrices and optical thickness in the range 0.25 to 64, it was concluded 
that ''approximations, based on polarization by single scattering only ... lead to large errors. Adding 
polarization by second order scattering gives only a small improvement" for the intensity and degree of 
linear polarization {Hovenier, 1971), 

Kawabata and Ueno (1988) developed algorithm for the vertically inhomogeneous media to 
compute the scalar intensity using the first three orders of scattering. This algorithm was advanced by 
Natraj and Spurr (2007) who provided the vector solution for the first two orders of scattering using the 
generalized spherical functions formalism (Gelfand et al., 1963; Siewert, 2000). Natraj and Spurr (2007) 
gave a detailed accuracy study of the approximate solution in application to the oxygen A-band, reporting 
that "results are exact in the center of strong lines and most inaccurate (-30% error in the Stokes 
parameter Q) in the continuum" . 

Although the theoretical derivations of the low orders of scattering VRTE approximation were 
previously reported, their accuracy evaluation for the aerosol atmospheres in the atmospheric window 
channels has been missing. This paper presents an accuracy study for an arbitrary anisotropic scattering 
matrix and both reflected and transmitted radiation. We pursue a practical target of 0. 2-0.5% in 
computation of polarized components of intensity. This level of accuracy dramatically limits the possible 
range of optical depths to Dave's (1964) assessment of 0.05. For this reason, this paper only considers the 
homogeneous case although the computational method could be generalized to include atmospheric 
vertical structure. For the same reason, the main accuracy analysis focuses on a conservative scattering 
case which provides the strongest constraint on the range of optical depths given the target accuracy. 

This paper is structured as follows: the statement of the problem, the algorithm of solution and 
computational relations, used in the paper, are given in Section (2). The definition of the cases, numerical 
results and discussion are presented in Section (3). The paper is concluded with a brief summary. 

2. Computational Relations 

We are considering the VRTE boundary problem without thermal sources for the homogeneous 
plane-parallel atmosphere over a black surface (Chandrasekhar, 1950). The VRTE solution is 
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decomposed into azimuthal Fourier series following Siewert (2000). The single scattering (SS) solution 
for the m-th harmonic can be obtained following van de Hulst (1948): 
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The following notations are used above; L = [/ Q U V] is the Stokes vector, L,, is the initial 

vector for solar irradiance, the overhead symbol indicates square matrices, x (0 < x < tq) is the 
optical thickness, and coo is the single scattering albedo (SSA). The directions are given by the solar 
(SZA, 0o) and view (VZA, 0) zenith angles (p=cos0). The positive (p^ = p>0) and negative directions 
correspond to transmitted and reflected light, respectively. The m-th Fourier term of the phase matrix is 

P” ((1, n') = y; (2«: + Dfi” (p)ft fl" (n') , (2) 

k=m 

where is defined by the greek-constants {de Rooij and van der Stap, 1984). n^(p) contains F(”(p) , 


/?™(p) and T/’'(p) polynomials related to generalized spherical functions (Gelfand et ah, 1963; Siewert, 
2000). K is the maximal considered order of the polynomials. The azimuth-dependent solution is obtained 
from Eq.(l) after summation Fourier series over m. 

Using Eq.(l) as a source function in VRTE and applying analytical integration over x (Hovenier, 

1971), the following result for the second order scattering L 2 (x,p) can be obtained: 

• Reflected radiance: 
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• Transmitted radiance: 
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In the developed code, a gaussian quadrature is used to perform integration in Eqs. (3) and (4). The 
singular points 1/ (p.'-p,)for p. — >p^ can be handled either by changing the order of quadrature to avoid 
view directions or by expanding the exponents in Taylor series. Below, the term “double scattered” (DS) 
will be used to describe contribution from the first two orders of scattering, L7(x,)Li)-f L 2 (x,)li) (Eqs. (1), 
(3), (4)). 


3. Numerical Results and Discussion 

The accuracy of SS and DS solutions was investigated using three types of scattering matrices 
with different scattering anisotropy. The first one is Rayleigh scattering matrix {Coidson et al, 1960) 
which has only two non zero expansion moments. The second one represents the fine spherical non- 
absorbing aerosol particles described by the gamma distribution with effective radius 0.2 pm, variance 
0.07 pm, and real refractive index 1.44 {de Rooij and van der Stap, 1984). The scattering matrix at the 
wavelength 0.951 pm has the average scattering cosine ASC = [p*.=i]j j =0.485 (see Eq.(2)) and A' = 11 

expansion terms. It was used by Evans and Stephens (1991) and is provided at http://nit.colorado.edu/. 
The third one, used in the code comparison study of Kokhanovsky et al. (2010) and available at 
www.iup.physik.uni-bremen.de/~alexk, represents the coarse aerosol fraction described by a lognormal 
size distribution with effective radius 0.3 pm, variance 0.92 pm, and real refractive index 1.385. It is 
characterized by the ASC = 0.793 and A' = 918 expansion moments. All K moments were kept in our 
computations. The phase matrices for aerosol cases are shown in Eigure 1. 

The accuracy of the result can be estimated using the relative error, e.g. =(I -l)\00/l [%], 

where / and / are the accurate and the approximate solutions, respectively. Because the polarization 
components Q, U, V may take zero values at certain points whose location is a priori unknown, we will 
further use the Chandrasekhar's vector with the first two components defined as // = (/ -f Q)!2, A = (/ - 
Q)H. Ii and A represent the parallel and perpendicular components of the polarized light intensity 
{Chandrasekhar, 1950; van de Hulst, 1957) which are non-negative. A may take a zero value only in the 
case of 100% linear polarization that rarely happens in the atmosphere. Thus, the accuracy is evaluated 
below using the following relative error [%]: 

'7)100//,,.. (5) 
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The accuracy of SS and DS solutions in this work was established by comparison with the vector 
code MVDOM {Budak et al, 2010). The relative errors Ei and Er for the case of conservative Rayleigh 
scattering for the reflected and transmitted radiance were studied as a function of VZA for two solar 
zenith angles SZA = 20°, 60° and five relative azimuths (p = 0°, 45°, 90°, 135°, 180°. Two optical depths x 
= 0.01 and x = 0.05 were chosen as in Dave (1964). The results obtained are in good agreement with those 
from the previous studies (Dave, 1964). At low optical depth, the error of SS approximation can be as 
high as several percent. It reaches 100% at 90° to the light incidence direction for the parallel component 
of intensity. In this direction, the amplitude of scattering of 7/ component is zero in the first order of 
scattering, and thus all of 7/ signal is created by the multiple scattering. This can be shown formally using 
the scattering matrix in Chandrasekhar's basis, which is given by the following expression 
{Chandrasekhar, 1950, Chapter 1, Eq.(201)) 

Pc,{y) = ^-^diag\_f 1 y y]- (6) 

At 90° scattering angle, its cosine is zero (y = 0) which cancels the first order scattering contribution for 
the 7/ component. With addition of the second scattering order, the error drops below 1-1.5% except when 
close to the 90° scattering angle when the maximal error reaches 12-18% at x=0.01 and 0.05, respectively. 

The error analysis for the fine aerosol fraction is presented in Figure 2 for the same view 
geometry. The top and bottom plots show the reflected and transmitted radiation, respectively. The dash 
and solid lines correspond to SZA=20° and 60°, respectively, and different azimuths are indicated by 
different colors. Compared to the Rayleigh case, the errors are a factor of ~2 higher at the same optical 
thickness. The average errors Er and Ei of the DS solution are approximately 0.6-0. 8% at x=0.01 and 1.8- 
2.5% at x=0.05, while the errors of the SS solution are about a factor of 8 higher. The maximal error for 
the parallel component of light is shifted from 90° towards larger scattering angles, which follows from 
the Mie theory {van de Hulst, 1957, p. 147, figure 23). 

Figure 3 repeats the error analysis for the coarse aerosol model. In this case, N = 600 ordinates 
were used in computations with = 918 expansion terms. The error decreases as compared to both the 
fine aerosol model (Figure 2) and Rayleigh case. It still remains at the level of several percent in the SS 
approximation, but drops down to 0.2-0. 3% and 1-1.5% in the DS approximation at x=0.01 and 0.05, 
respectively. 

The above analysis was conducted for the conservative scattering case. An aerosol absoiption 
reduces the number of scattering events thereby increasing the accuracy of SS and DS approximations. 
However, the error changes insignificantly within the realistic range of aerosol SSA values 0.8-1. 

The analyzed approximate solution does not describe well the neutral points of polarization of the 
clear sky {Chandrasekhar, 1950). These points appear only in the DS approximation, but are shifted with 
respect to the true position {Hovenier, 1971). For example, at x=0.2, the Babinet and Brewster points (cp = 
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0“) are shifted towards incidence direction by approximately 3“ and 5°, respectively, while the Arago 
point (SZA=78.5° and cp=180°) is shifted 5° towards the horizon. 

The results presented above allow us to extend conclusion of Dave, made for the Rayleigh case, to 
the aerosol atmospheres. Namely, the first two scattering orders provide the high accuracy of polarization 
computations ~0. 2-0.5% only when the media is optically thin (t< 0.01) or there is a significant gaseous 
absorption which essentially eliminates higher orders of scattering (Natraj and Spun, 2007). In the 
typical case of aerosol remote sensing in the atmospheric window channels, the accuracy of the DS 
solution deteriorates rapidly with the optical thickness making it impractical for the data analysis. 

One interesting finding of this analysis is a weak dependence of the relative errors Ei and Er, only 
within a factor of 2-4, among the three considered cases with the large difference in the scattering matrix 
anisotropy. One explanation is that all considered examples were generated for the optically thin media 
(x<0.1) where the first two scattering orders dominate the signal regardless of the scattering anisotropy. 
To confirm this statement, a numerical experiment was performed with the Monte-Carlo code (Marshak 
and Davis, p. 264) simulating the transfer of photons along a single line (ID transport). The photon can 
scatter forward with the probability p+=0.5-(l+g) or backward with the probability p_=0.5-(l - g), where g 
is the asymmetry parameter, or ASC. Figure 4 shows the number of photons scattered once (n\, blue line) 
and twice (ni, green line) and their sum (nj -f U 2 , red line) as a fraction of the total number of scattering 
events (N). These fractions correspond to the reflected (transmitted) radiation calculated using the first 
order or the two orders of scattering. The superscripts ‘V” and “r” denote the reflected and transmitted 
radiation, which was computed for SSA=1 and three values of optical depth i = 0.1, 1, 10. The ratio 
between the forward and backward scattering (1 -i- g) / (1 - g) ranges from 5.67 to 9, 19, 199 for g = 0.7, 
0.8, 0.9, 0.99 respectively. The left plots of Figure 4 at x=0.1 shows that the relative contribution of the 
first two scattering orders (ni/N and no/N) changes little while the asymmetry of phase function changes 
significantly with g. This, however, is not the case when higher orders of scattering become prominent at 
larger optical depths x=l, 10. Note that for larger optical depths (right panels), the single scattering 
contribution to the total reflected or transmitted radiation decreases with the increase of the phase 
function anisotropy. 

All of our previous results considered an optically thin case. At the opposite end of very high 
optical depths, the accuracy of SS and DS solutions for Q and U components rises again because the high 
scattering orders with alternating sign from rotation of the reference plane effectively cancel each other. 
This effect is illustrated in Figure 5 showing Stokes components for radiation reflected from a cloud with 
x=50 (SSA=1). The “cloud” model (Kokhanovsky et al., 2010), available at www.iup.phvsik.uni- 
bremen.de/~alexk , is characterized by an effective radius 5 pm, variance 0.4 pm, and maximal droplet 
size of 100 pm. Although the approximate result has a rather low accuracy, the DS-solution captures a 
large part of Q and U components, unlike / and V which require accumulation of all orders of scattering. 
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4. Conclusion 

This paper studied the accuracy for the first two scattering orders using code MVDOM as a 
reference. Although the previous studies conducted for Rayleigh scattering {Dave, 1963) showed that the 
DS approximation has a good accuracy only at low optical depths (t < 0.05), it was not clear whether this 
conclusion holds for the aerosol atmospheres with higher scattering anisotropy which tends to reduce the 
effect of light polarization on the total intensity (Mishchenko et al., 2006). Our analysis, which included 
two aerosol scattering matrices with weak and strong anisotropy, corroborated Dave's conclusion showing 
that DS approximation provides an accuracy of better than 0.2-0.5% at t~ 0.01, and 1-2% at t~ 0.05, 
however a full multiple scattering solution is required at higher optical depths. At low optical depths, this 
error has a weak dependence on the anisotropy of the scattering matrix. 

Although we demonstrated that the DS approximation cannot be used at t > 0.05, it has a certain 
didactic value as an approximate solution for the degree of polarization at the near 90° scattering angles 
as well as for the V-component of the Stokes vector. In these cases, the primary scattering gives a zero 
contribution (e.g., see Figure 5), and the signal is created by the higher orders of scattering. Presently, 
both theory and passive remote sensing applications pay little attention to the atmospheric ellipticity as V- 
signal is much weaker than the other components of the Stokes vector. On the other hand, as well as the 
total intensity, F-component does not depend on the choice of the reference plane, and its information 
content is not reduced by the higher orders of scattering contrary to the Q and (7 components. The DS 
approximation thus may serve as a useful constraint for the V-component in accurate computations. 
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Figure 1. Phase matrices for aerosol cases defined in the paper. The maximal values of the phase 
functions ai of the fine and coarse fractions are approximately 4 and 1500 respectively; aj, as, 
hi, h 2 are the [1, 1], [3, 3], [1, 2] and [3, 4] elements of the phase matrix respectively. 

Figure 2. Relative error (%) for the parallel (Figures 2(a) and 2(b)) and perpendicular (Figures 
2(c) and 2(d)) components of the reflected (VZA = 90°...180°) and transmitted (VZA = 0°...90°) 
radiation for the fine aerosol model. The error is computed for five different azimuthal angles 
(p=0“, 45°, 90°, 135°, 180° shown by different colors (refer to the legend, Figure 2(e)), and two 
solar zenith angles SZA=20" (dash lines) and 60° (solid lines). The results are shown for two 
values of optical thickness, t = 0.01 and t = 0.05 at SSA=1. Figures 2(a) and 2(c) show the result 
for the single scattering approximation. Figures 2(b) and 2(d) show the result for the double scat- 
tering approximation. 

Figure 3. The same as in Figure 2 but for the coarse aerosol fraction. The legend is shown in Fu- 
gure 2(e). 

Figure 4. The Monte-Carlo simulations of the number of single scattered photons (ni, blue line), 
re-scattered photons (ri 2 , green line) and their sum («/ -i- ri 2 , red line) as a fraction of all scattered 
photons (A^). SSA = 1 (see sec, 3). The superscripts “r” and “r” refer to the reflected and trans- 
mitted photons respectively. The ASC is the average scattering cosine. 

Figure 5. 1, Q, U, and V components of the Stokes vector of reflected light for a cloud {Kokha- 
novsky et ah, 2010) with x = 50 (see sec. 3). The view geometry is described by SZA=30°, (p=45° 
(Fig. 5a), and cp=135° (Fig. 5b). The green, blue and red lines represent the single scattering (SS), 
double scattering (DS), and the full multiple scattering (MS) solutions. 
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1. The double scattering model has accuracy below 0.2-0.5% at T < 0.01 in the conservative case 

2. The error rapidly grows with optical depth, so the full VRTE solution is required for the 
practical aerosol remote sensing; 

3. At low optical depths, the relative error for low scattering orders weakly depends on the 
scattering anisotropy. 




